% This code computes initial rate regulation costs using the FOC w.r.t. p1.
% cost1_e: cost coefficient, NN x 1
% cost_ml: total min loss regulation cost, NN x 1

function [cost1_e, cost_ml] = minloss_estimation(lr_target, cost2_1, ...
    D_input, ...
    resource_y1, resource_y2, pr_y, N_y, ...
    K, M,...
    myopic, delta, B1_c, B2_c, B1_f, B2_f, ...
    beta_c, beta_f, T1, alpha, rho, crra, sigma_c1, Er, cost_weight, use_cost_weight)
%% Get data
s1_predicted = D_input.s1_predicted;
s1_predicted_y_c = D_input.s1_predicted_y_c;
s1_predicted_y_m = D_input.s1_predicted_y_m;
p1 = D_input.p1;
p2 = D_input.p2;
pr_s = D_input.pr_s;
claims2 = D_input.claims2;
delta_p = D_input.delta_p;
major = D_input.major;
mkt = D_input.mkt;
in_regression = D_input.in_regression;

%% Calibrate p1_target (NN x 1). 
p1_target = p1_target_ftn(lr_target,  pr_s, claims2, ...
    beta_f, T1, B2_f,  B1_f, Er);

%% Market share derivatives w.r.t. p1
NN = size(cost2_1,1);

[D_s1_p1, D_s2_p1, ~, D_s1_ind_p1]...
    = D_s_price(p1, s1_predicted_y_c, s1_predicted_y_m, pr_s, ...
    resource_y1, resource_y2, pr_y, N_y, ...
    NN, K, M,...
    myopic, delta, B1_c, B2_c, beta_c, T1, alpha, rho, crra, sigma_c1);

%% Compute reg costs using FOC w.r.t. p1
D_profit1_p1 = B1_f*(s1_predicted + p1.*D_s1_p1);

KM = K*M;
D_claim_av = zeros(NN, KM);
pr_stay_fixed2 = (1-unique(delta))*ones(N_y,KM);
for j=1:NN
    D_s1_ind_p1_jj = D_s1_ind_p1(:,j); % N_y x 1
    D_s2_ind_p1_jj = repmat(D_s1_ind_p1_jj,1,KM).* pr_stay_fixed2; %(N_y, KM) matrix
    for k=1:KM
        D_claims2_jj = claims2(j,k); % scalar
        D_claim_av(j,k) = sum(cost_weight.*repmat(D_claims2_jj, 1, N_y)'.*D_s2_ind_p1_jj(:,k).*pr_y, 'all');
    end
end

if use_cost_weight==0
    D_profit2_p1_k = (p2 - claims2).*D_s2_p1 + (delta_p>0).*cost2_1.*delta_p;
else
    D_profit2_p1_k = p2.*D_s2_p1 - D_claim_av + (delta_p>0).*cost2_1.*delta_p;
end

D_profit2_p1 = (beta_f^T1)*B2_f*sum(D_profit2_p1_k.*pr_s,2);

% profit derivative w.r.t. p1, N.ist x 1 
D_profit = D_profit1_p1 + D_profit2_p1;

% deviation from p1_target, N.ist x 1 
p1_deviation = p1-p1_target;

% coefficient on the 1st-period min ratio regulation cost, NN x 1 
cost1_e = (1./p1_deviation).*D_profit;

% total 1st-period min ratio regulation cost, NN x 1 
cost_ml = 0.5*cost1_e.*(p1_deviation.^2);
